model{
  for (i in 1:n) {
    for (j in 1:m) {
      for (k in 1:t) {
        y[i, j, k] ~ dnorm(mu[i, j, k], inv.phi2[j])
        mu[i, j, k] <- beta.raw[j, k] * theta.raw[i, k]
      }
    }
  }
  for (i in 1:n) {
    theta.raw[i, 1] ~ dnorm(0, inv.xi2)
    for (k in 2:t) {
      theta.raw[i, k] ~ dnorm(theta.raw[i, k - 1], inv.tau2.raw[k - 1])
    }
    for (k in 1:t) {
      theta[i, k] <- (step(beta.raw[1, k]) - (1 - step(beta.raw[1, k]))) *
        (theta.raw[i, k] * inv.xi)
    }
  }
  for (j in 1:m) {
    beta.raw[j, 1] ~ dnorm(0, 1)
    for (k in 2:t) {
      beta.raw[j, k] ~ dnorm(beta.raw[j, k - 1], inv.omega2.raw[k - 1])
    }
    for (k in 1:t) {
      beta[j, k] <- (step(beta.raw[1, k]) - (1 - step(beta.raw[1, k]))) *
        (beta.raw[j, k] * xi)
    }
    inv.phi2[j] ~ dgamma(1, 0.2)
  }
  phi <- 1 / sqrt(inv.phi2)
  inv.xi2 ~ dgamma(0.5, 0.5)
  inv.xi <- sqrt(inv.xi2)
  xi <- 1 / inv.xi
  kappa[1] ~ dunif(0, 1)
  kappa[2] ~ dunif(0, 10)
  lambda[1] ~ dunif(0, 1)
  lambda[2] ~ dunif(0, 10)
  for (k in usual.period) {
    inv.tau2.raw[k] <- inv.xi2 / (kappa[1] * kappa[1])
    inv.omega2.raw[k] <- 1 / (lambda[1] * lambda[1])
  }
  for (k in unusual.period) {
    inv.tau2.raw[k] <- inv.xi2 / (kappa[1] * kappa[1] * kappa[2] * kappa[2])
    inv.omega2.raw[k] <- 1 / (lambda[1] * lambda[1] * lambda[2] * lambda[2])
  }
  tau[1] <- kappa[1]
  tau[2] <- kappa[1] * kappa[2]
  omega[1] <- lambda[1] * xi
  omega[2] <- lambda[1] * lambda[2] * xi
}